rm(list=ls())
library(foreign)
library(locfit)
library(Hmisc)
#setwd("/Users/innoshuadmin/Box Sync/DRD")
setwd("/home/users/innoshu.AD3/drdfinal")
source("./functions/DRDfunctions.R")
## load in robust bandwidth selection functions from CCT(2014)## 
source("./functions/rdbwselect.R")
source("./functions/functions.R")

## a list of functions used in the following codes are
## 1. rdbwselect -- bandwidth selection function, from CCT (2014)
## 2. MRD_sharp -- mean RD regression, sharp design or reduced form MRD regresion only. 
## 3. DRD_test -- conducts distribution RD test proposed in this paper, calls for function DRD_sharp.
####### 3.1. DRD_sharp -- distribution RD regression, sharp design or reduced form MRD regresion only. 
## 4. MRD_plot -- generates mean RD plots used in this paper, uses results from MRD_sharp.
## 5. DRD_plot -- generates distributional RD plots used in this paper, uses results from DRD_sharp.
##                Two types pf DRD figures can be generated. type=1 produces reduced form DRD estimates with pointwise CI;
##                type=2 produces scaled DRD estimates and compare them with uniform critical values.
## 6. DRDpvalue -- tabulates p-values.

## load in data
romanian <- read.dta(file="./APP1/romanian_2sch.dta")
newdata<- romanian[c("bcg_Dmean","agus","dzag","zga","max_z")]
newdata <-subset(newdata,is.na(newdata$bcg_Dmean)==0&dzag!=0)
newdata <-subset(newdata,abs(dzag)<=min(max(dzag),max(-dzag)))

## following are the main codes
numd=3
pcresults=matrix(NA,ncol=5,nrow=numd)
mrdresults=matrix(NA,ncol=5,nrow=numd)
drdresults=matrix(NA,ncol=7,nrow=numd)
ctff<-0
Xeval <- seq(-0.5, 0.5, length=51)  

for (table in 1:3){
print(paste("Table ",table,sep=""))
for (d in 1:numd){
	if (d==1) {data<-newdata} else
	if (d==2) {data<-subset(newdata,zga>=7.37)} else
	if (d==3) {data<-subset(newdata,zga<7.37)} 

#################################################################
daty=data$bcg_Dmean
datx=data$dzag
datt=data$agus
print(summary(datx))
N=length(daty)
ygrid=quantile(daty,seq(0.001,0.999,length.out=999))

if (table==1){
### the lines below are for reproducing Table 1
if (d==1){
smooth=4.5
mrdbw.out = rdbwselect(daty,datx, bwselect = "CCT")
bw_mrd=mrdbw.out$bws[1,1]*N^(1/5-1/smooth)
}
## these lines above are for reproducing Table 1
} else
if (table==2){
### the lines below are for reproducing Table 12
smooth=4.5
mrdbw.out = rdbwselect(daty,datx, bwselect = "CCT")
bw_mrd=mrdbw.out$bws[1,1]*N^(1/5-1/smooth)
## these lines above are for reproducing Table 12
} else
if (table==3){
### the lines below are for reproducing Table 13
if (d==1){
smooth=4
mrdbw.out = rdbwselect(daty,datx, bwselect = "IK")
bw_mrd=mrdbw.out$bws[1,1]*N^(1/5-1/smooth)
}
## these lines above are for reproducing Table 13
} else
if (table==4){
### the lines below are for a robustness check table that is not reported in the paper
smooth=4
mrdbw.out = rdbwselect(daty,datx, bwselect = "IK")
bw_mrd=mrdbw.out$bws[1,1]*N^(1/5-1/smooth)
## these lines above are for a robustness check table that is not reported in the paper
}

## proportion of complier
results<-MRD_sharp(datt,datx,bw_mrd,Xeval=Xeval)
pcresults[d,1:5]<-round(c(bw_mrd,results$estimate,(results$estimate/results$tstat),(2*pnorm(-abs(results$tstat))),sum(abs(datx)<=bw_mrd)),3)

## Mean Change in Outcome
results<-MRD_sharp(daty,datx,bw_mrd,Xeval)
mrdresults[d,1:5]<-round(c(bw_mrd,results$estimate,(results$estimate/results$tstat),(2*pnorm(-abs(results$tstat))),sum(abs(datx)<=bw_mrd)),3)

if (table==1){
# MRD plot
xlab="Standardized Running Variable"
ylab1="Baccalaureate Grade"
ylab2="# of Individuals"
main="Mean RD / Histogram"
labels<-data.frame(xlab,ylab1,ylab2,main)
results$bw=bw_mrd
pdf(paste("./APP1/figures/romanian_Figure_MRD_CCT_d",d,".pdf",sep=""),width=7,height=5)
MRDplot(Xeval,daty,datx,results,labels)
dev.off()
}

## DRD test
if (table<=2){
results<-DRDtest(daty,datx,ygrid, h=bw_mrd, bwsel = "CCT",undersmooth=smooth)
} else
if (table==3){
results<-DRDtest(daty,datx,ygrid, h=bw_mrd, bwsel = "IK",undersmooth=smooth)
}
drdresults[d,1:7]<-round(c(results$bw,results$bw_iter,min(results$ksstat_iter),max(results$ksstat_iter),DRDpvalue(min(results$ksstat_iter),side=1),DRDpvalue(max(results$ksstat_iter),side=1),sum(abs(datx)<=results$bw_iter)),3)

if (table==1){
# DRD plot
xlab="Baccalaureate Grade"
ylab1="Conditional CDFs"
ylab2="Difference in CCDFs"
main="Distributional RD"

labels<-data.frame(xlab,ylab1,ylab2,main)
pdf(paste("./APP1/figures/romanian_Figure_DRD_CCT_d",d,".pdf",sep=""),width=7,height=5)
DRDplot(ygrid,results,labels,iter=1)
dev.off()
}
}
print("First stage: bw-estimate-sterr-2sidedpvalue -- nobs")
print(pcresults)
print("Reduced form MRD: bw-estimate-sterr-2sidedpvalue-stat -- nobs")
print(mrdresults)
print("Reduced form DRD: bw-iterated bw-statmin-statmax-statminpvalue-statmaxpvalue-- nobs")
print(drdresults)
}
